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By J. E. Taylor^ and K. J. Worsley^ 

Stanford University, Universite de Montreal and McGill University 

Our data are random fields of multivariate Gaussian observations, 
and we fit a multivariate linear model with common design matrix at 
each point. We are interested in detecting those points where some 
of the coefficients are nonzero using classical multivariate statistics 
evaluated at each point. The problem is to find the P-value of the 
maximum of such a random field of test statistics. We approximate 
this by the expected Euler characteristic of the excursion set. Our 
main result is a very simple method for calculating this, which not 
only gives us the previous result of Cao and Worsley [Ann. Statist. 
27 (1999) 925-942] for Hotelling's T^, but also random fields of Roy's 
maximum root, maximum canonical correlations [Ann. Appl. Probab. 
9 (1999) 1021-1057], multilinear forms [Ann. Statist. 29 (2001) 328- 
371], [Statist. Probab. Lett 32 (1997) 367-376, Ann. Statist. 25 
(1997) 2368-2387] and scale space [Adv. m Appl. Probab. 33 (2001) 
773-793]. The trick involves approaching the problem from the point 
of view of Roy's union-intersection principle. The results are applied 
to a problem in shape analysis where we look for brain damage due 
to nonmissilc trauma. 

1. Introduction. Our motivation comes from a study by Tomaiuolo, Wors- 
ley, Lerch, Di Paulo, Carlesimo, Bonanni, Caltagirone and Paus [25] on the 
anatomy of a group of 17 nonmissile brain trauma patients measured by 
magnetic resonance imaging (MRI). The aim is to detect regions of brain 
damage (shape change) by comparing anatomy at each point in = 3 di- 
mensional space to that of a group of 19 age and sex matched controls. Each 
brain was first linearly transformed into a common stereotactical reference 
space. Then the method of Collins, Holmes, Peters and Evans [6] was used 
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to find the nonlinear vector deformations Yi{s) G JR'^ (d = 3 here) required to 
transform the MRI image of subject i to a common atlas standard at each 
point s inside a search region S C usually the whole brain — see Figure 
1. This sort of anatomical data is good at detecting changes in the bound- 
ary of brain structures, with the added benefit of estimating the direction 
in which the change took place. 

To do this, Cao and Worsley [4] set up a linear model for subject i, 
i = 1, . . . ,n: 

(1) nsy = x',(3{s) + Zi{syj:{s)'/^, 

where Xi is a p- vector of known regressors, and /?(s) is an unknown p x d 
coefficient matrix. The error Zi{s) is a d- vector of independent zero mean, 
unit variance Gaussian random field components with the same spatial cor- 
relation structure, and Var(yi(s)) = S(s) is an unknown d x d matrix. We 
can now detect how the regressors are related to shape at point s by testing 
contrasts in the rows of f3{s). Classical multivariate test statistics evaluated 
at each point s then form a random field T(s). 

We suspect that the changes are confined to a small number of localized 
regions and the rest of the brain is unchanged. It is not hard to show that 
if the spatial pattern of the change matches the spatial correlation function 
of the errors, and if T(s) is the likelihood ratio at a single point s, then the 
spatial maximum of T{s) is the likelihood ratio test under unknown signal 
location, provided is known (Siegmund and Worsley [14]). In essence, 
this is just the Matched Filter Theorem from signal processing. Thresh- 
olding T{s) at that threshold which controls the P-value of its maximum 
should then be powerful at detecting changes, while controlling the false 
positive rate outside the changed region to something slightly smaller than 
the nominal P-value. Our main problem is therefore to find the P-value of 
the maximum of random fields of multivariate test statistics. 

The outline of the paper is as follows. An approximate P- value of the 
maximum of a random field is given in Section 2, and in Section 3 we apply 
this to random fields of multivariate test statistics T{s) such as Hotelling's 
T^, Roy's maximum root and maximum canonical correlation. The same 
methods can also be applied to random fields of (Lin and Lindsay [12] 
and Takemura and Kuriki [18]) statistics, but this will be the subject of a 
forthcoming paper (Taylor and Worsley [23]). In Section 4 we apply these 
results to the nonmissile brain trauma data above. A further application to 
scale space random fields is given in Appendix A. 5. 
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2. P-value of the maximum of a random field. 



2.1. The expected Euler characteristic. A random field T(s), s G 5" C M^, 
is isotropic if it has the same distribution as T{a + Ss), where a E is 
any translation and B is any orthnormal (rotation) matrix. A very accurate 
approximation to the P-value of the maximum of such a random field, at high 
thresholds t, is the expected Euler characteristic (EC) ip of the excursion set 

(2) pLaxr(s)> A ^¥.{^{seS:T{s)>t]) = Y,^H{S)pi{t), 

where fJ.i{S) is the i-dimensional intrinsic volume of S, and Pi{t) is the 
i-dimensional EC density of the random field above t (Adler [1, 2], Wors- 
ley [27], Taylor, Takemura and Adler [22] and Adler and Taylor [3]). The 
accuracy of the approximation (2) will be discussed in Section 2.2. 
For = 3, our main interest in applications, the EC of a set is 

if = tjblobs — jjhandles or tunnels + (jinterior hollows 

(see Figure 2). It can be evaluated numerically for a subset of a rectilinear 
mesh by 

(3) if = ttpoints — Jjedges + jjfaces — ftcubes, 

where, for example, a cube is a set of 8 adjacent mesh points, differing by one 
mesh step along each axis, and all inside the set [Figure 2(a)]. This method 
was used to calculate the EC of the excursion sets in Figure 2(b)-(e). 
Intrinsic volumes are defined in Appendix A.l, and for = 3, they are 

Ho{S) = (p{S), 

/ii(5) = 2 caliper diameter (S), 

(4) 

A*2(5') = 2 surface area(5), 
IJ-siS) = volume(5). 

For convex S, the caliper diameter is the distance between two parallel 
planes tangent to S, averaged over all rotations. 

EC density is defined in Section 2.3. As an example, suppose T = T{s) is 
a unit Gaussian random field (UGRF), defined as a Gaussian random field 
with E(r) = 0, Var(r) = 1 and Var(r) = InxN, the N x N identity matrix, 
where dot denotes (spatial) derivative with respect to s. Then the first four 
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(0 ftij 

Fig. 1. Shape analysis of nonmissile brain trauma data, (a) Trauma minus control av- 
erage deformations (arrows and color bar), sampled every 6mm inside the brain, with 
T{s) = Hotelling' s field for significant group differences (threshold t = 54.0, P = 0.05 ). 
The reference point of maximum Hotelling' s T is marked by the intersection of the 
three axes, (b) Closeup of (a) showing that the damage is an outward movement of the 
anatomy, either due to swelling of the ventricles or atrophy of the surrounding white mat- 
ter, (c) Regions of effective anatomical connectivity with the reference point, assessed by 
T{s) = maximum canonical correlation field (threshold t = 0.746, P = 0.05j. The reference 
point is connected with its neighbors (due to smoothness) and with contralateral regions 
(due to symmetry), (d) Regions where the connectivity is different between trauma and con- 
trol groups, assessed by T{s) = Roy's maximum root field (threshold t — 30.3, P = 0.05 ). 
The small region in the contralateral hemisphere (right) is more correlated in the trauma 
group than the control group. 



EC densities of T are 

Poit)=¥{T>t), 
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(a) (d) (e) 

Fig. 2. (a) The Euler characteristic of a 3D set is 
(p — ttfeZofes — fihandles or tunnels + '^hollows = 2 — 1 + = 1, or, if we fill the set with a 
fine rectilinear mesh, then ip = '^points ~ '^edges + 'if aces — ^cubes = 47 — 70 + 26 — 2 = 1. 
(b) The excursion set of a 3D isotropic Gaussian random field with mean zero and 
variance one above a threshold t = —2; the set contains isolated hollows that each 
contribute +1 to give tp = 5 with E(</3) = 6.7 from (2); (c) at t = the handles or tunnels 
dominate, each contributing —1 to give (p — —28 with E((^) = —20.0; (d) at t — 2 the 
handles and hollows disappear, leaving isolated blobs, each contributing +1 to give (p — 11 
with E(<^) = 16.1; (e) at t — 3 only one blob remains (containing the maximum value of 
3.16) to give <p = 1 with E(</5) = 2.1. At very high thresholds K{<p) is a good approximation 
to the P-value of the maximum. 



pi(t) = (27r)-iexp(-tV2), 

(5) 

/j2(t) = (27r)-3/2texp(-tV2), 

P3(t) = (27r)-2(t2_i)exp(-tV2). 

Note that any stationary Gaussian random field can be transformed to a 
UGRF by an appropriate linear transformation of its domain and range. In 
particular, if Var(T) = cInxN for some scalar c, then pi{t) is multiplied by 



2.2. The accuracy of the approximation. A heuristic explanation for why 
we use the expected EC as a P-value approximation is as follows. If the 
threshold t is high, the handles and hollows of the excursion set tend to dis- 
appear, leaving a set of isolated blobs, each containing one local maximum, 
so that the EC then counts the number of connected components [Figure 
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2(d)]. At very high thresholds the excursion set is mostly empty with an EC 
of 0, or occasionally, when maxs^s T{s) > t, it will contain just one connected 
component with an EC of 1 [Figure 2(e)]. Thus, at these high thresholds the 
expected EC is a good approximation to the P- value of maxg^s T (s) . The 
beauty of the EC is that there is an exact expression for its expectation for 
all thresholds. 

Moreover, the approximation (2) is astonishingly accurate when S is either 
convex or has smooth boundary, and T{s) is Gaussian, in which case 

E{ip{s£S:T{s)>t}) 

(6) 

= P(r >t) + {co + cit + --- + CD-it^"^) exp(-tV2) 

for some constants cq, . . . ,cd-i [see (5) and (9) below]. Since P(T > t) = 
0(l/t) exp(— 1^/2), it might be thought that the error in the P-value ap- 
proximation (2) is simply the next term down, that is, 0(l/t^) exp(— 1^/2). 
In fact, the error is exponentially smaller than this: 

(7) p(^maxT(s) = E((^{s G 5 : T(s) > t}) + 0(exp(-atV2)) 

for some a> 1 which is related to the curvature of the boundary of S and 
Var(r) (Taylor, Takemura and Adler [22]). This means that there are in 
effect no further terms in the expansion, and the expected EC captures 
essentially all of the polynomial terms in the P-value. 

We will see later that many of the P-value approximations for non- 
Gaussian fields can be transformed into P-value approximations of Gaus- 
sian fields on larger parameter spaces. This means that these special non- 
Gaussian fields also have exponentially accurate P-value approximations as 
in (7). There will be some examples, however, with a field in the denom- 
inator, for which this trick will not work. For these examples, it is difficult 
to give quantitative bounds on the error. The general techniques in Taylor, 
Takemura and Adler [22] still apply to these non-Gaussian fields, though the 
final form of the bound is not explicitly known in these cases. 

2.3. The EC density. A direct method for finding EC densities of any 
sufficiently smooth random field can be derived using Morse theory (Morse 
and Cairns [13]) for i > 0, 

(8) piit) = E{l{T>t} det{-f,) 1 1, = 0)P(r, = 0), 

where dot notation with subscript i denotes differentiation with respect to 
the first i components of s, and double dot with subscript i denotes the 
matrix of second derivatives with respect to the first i components of s 
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(Adler [1] and Worsley [27]). For i = 0, po{t) = P(T > t). If T is a UGRF, 
Adler [1] uses this method to obtain its EC density 

(9) P?ii) = (^i\nT>t), 



which leads directly to (5) and (6). 

Our main interest in this paper is finding P-values for maxima of random 
fields of test statistics commonly encountered in multivariate analysis, with 
the ultimate aim of applying this to detecting points s where coefficients 
in a multivariate linear model are nonzero. First of all we must define the 
random field T. For example, a random field with d degrees of freedom 
is defined as 

(10) T{s) = Z{syZ{s), 

where Z{s) is a d-vector of i.i.d. UGRFs. In a similar way we can define 
other test statistic random fields, such as t and F statistic random fields, 
by applying the usual definition to component UGRFs. 

2.4. The Gaussian kinematic formula. Finding workable expressions for 
EC densities is a tedious process. The result for UGRFs (9) appears to be 
simple, but its derivation takes most of Chapter 5 of Adler [1]. In fact, Adler 
did not at first recognize that it could be expressed as the derivative of the 
P-value at a point. For many years this was thought to be a coincidence, 
but a new result of Taylor [20] on Gaussian kinematic formulas shows that 
indeed this is not a coincidence in certain cases (the Gaussian is one, the 
is another): the EC densities of functions of UGRFs fall out as the 
coefficients of a power series expansion of the probability content of a tube 
about the rejection region, in terms of the tube radius. For the Gaussian 
case, this tube is just another rejection region of the same shape, which 
ultimately leads to the special form of the Gaussian EC density (9). 

The details of the Gaussian kinematic formula are as follows. Suppose 
T(s) = f{Z{s)) is a function of Z{s), a d-vector of i.i.d. UGRFs as be- 
fore, so that the rejection or critical region is C = {z:f{z) >t}c W^. Let 
Tube(C, e) = {x:mmz£c \\^ — x\\ < e} be the tube of radius e around C. 
Let Z he a d-vector of i.i.d. N(0, 1) random variables. Then the Gaussian 
kinematic formula is 



(11) P(Z G Tube(C, £)) = J2 -(2vr)^/Vi(t). 

For T{s) a UGRF, d=l, f(z) = z, and this leads directly to the Gaussian EC 
density (9). For an F statistic random field with {ri,u) degrees of freedom, 

/(3)=4^, 
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where z = (zj, 2:2)') ^\ £ ^ ■, ^2 G and d = ij + v. The rejection region C is 
a cone, and after a httle elementary geometry, it can be seen that the first d 
terms come only from the sides of the cone and are easy to calculate (Taylor 
and Worsley [23]): 

(12) F{Z eTuhe{C,e)) =F{VV > ^Jwt^ - e^Jl + t^) + 0{£'^), 

where V ^ independently of ~ • Expanding the right-hand side of 
(12) in powers of e and equating it to that of the Gaussian kinematic formula 
(11) gives the i-dimensional EC densities (t) of the F field for i <d. This 
agrees with an explicit expression found using the Morse theory result (8) 
that can be found in Worsley [26]. 

The higher-order terms in 0{e'^) in (12) come from the apex of the cone. 
Luckily we do not have to calculate them, since we only need EC densities 
up to dimension A'^ and if A'^ > d, then the F field is not defined (Worsley 
[26]). The F field is not defined because both numerator and denominator of 
F can take the value (i.e., F = 0/0) with positive probability. The reason 
is that the zero set of any of the component UGRFs is a smooth surface of 
dimension — 1; the intersection of d such surfaces is a set of dimension 
N — d which is nonempty (with positive probability) \i N — d>Q. On this 
set both numerator and denominator of F take the value with positive 
probability and F is not defined. 

2.5. Non-Gaussian random fields. Unfortunately the EC densities of non- 
Gaussian random fields are not always derivatives of the P- value as in (9). 
Instead we must go back to the Morse theory result (8) or, if the random field 
is a function of i.i.d. UGRFs (as is the case here), we can use the Gaussian 
kinematic formula (11). Both are equally complex for the types of random 
fields we are interested in. 

For the Gaussian kinematic formula, we must evaluate the probability 
content of a tube about a very complex rejection region. For the Morse 
theory result, we must obtain expressions for the joint distribution of the 
random field and its first two derivatives, then evaluate the expectation 
(8) by careful manipulation of vector and matrix random variables. This 
has been done for , t and F random fields (Worsley [26]), Hotelling's 
random fields (Cao and Worsley [4]) and correlation random fields (Cao and 
Worsley [5]). In each case, after a great deal of tedious algebra, an exact 
closed-form expression for the EC density is obtained. 

The types of random fields we are interested in are generalizations of 
Hotelling's T^, which is used to test for a single coefficient in a multivariate 
linear model. Our goal is to generalize random field theory to testing for 
multiple coefficients, the only missing piece. The obvious choice is the likeli- 
hood ratio statistic, Wilks' A, but so far this has resisted our attempts. In- 
stead, by a very simple trick that builds on previous EC densities and avoids 
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any further evaluation of the Morse theory result (8), we can easily obtain 
P-value approximations similar to (2) (but not EC densities) for Roy's max- 
imum root, a common alternative to Wilks' A. Using the same method, we 
re-derive in Section 3.1 the EC density for Hotelling's T^, shortening the 
original derivation from an entire Annals paper (Cao and Worsley [4]) down 
to just several lines. We use the same trick to get the scale space field from 
the Gaussian scale space field, again reducing most of an Advances in Ap- 
plied Probability paper (Worsley [28]) down to just one line — see Appendix 
A.5. 

2.6. The union-intersection principle. The trick is to approach the prob- 
lem from the point of view of Roy's union-intersection principle. Take, for 
example, the random field with d degrees of freedom (10). We can write 
this as the maximum of linear combinations of d i.i.d. UGRFs Z{s) as fol- 
lows: 

T{s) = maxT{s, u), where T{s, u) = {Z{syu)'^, 

and G M'^ is restricted to the unit sphere Ud = {u: u'u = 1}. In other words, 
we have written the field in terms of the square of a Gaussian field over 
a larger domain S x Ud- 

At first glance this seems to make things more difficult, but, in fact, it 
makes things easier. We can now apply the P-value approximation (2) to 
T{s,u), replacing S hy S x Ud- Using integral geometry, the intrinsic volume 
of 5 X f/rf is a simple function of the intrinsic volumes of S and Ud- The EC 
densities of T{s,u) are easily obtained from the Gaussian EC densities (9). 
Putting them all together and pulling out the coefficient of Hi{S) gives us 
the EC density for the field. 

2.7. Nonisotropic fields and Lipschitz-Killing curvature. There is only 
one technical difficulty. Although T{s,u) is isotropic in s for fixed u, it is 
not isotropic in u for fixed s (since u is on a sphere), nor jointly isotropic 
in {s,u). To handle this situation, Taylor and Adler [21] and Taylor [20] 
have extended the theory of Adler [1] and Worsley [27] to nonisotropic fields 
on smooth manifolds with piece- wise smooth boundaries. They show that 
if the random field is a function of i.i.d. nonisotropic UGRFs, then it is 
only necessary to replace intrinsic volume /Uj(S') in the expected EC (2) by 
Lipschitz-Killing curvature (LKC) Ci{S): 

(13) P (maxT{s) >t]^ K{(p{s G S : T{s) >t})=Y^ £i{S)pi{t). 

This result holds for nearly all nonisotropic random fields that are N(0, 1) 
at each point and smooth, that is, with at least two almost sure derivatives 
and an additional mild regularity condition; see Taylor and Adler [21]. 
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The LKC of is a measure of the intrinsic volume of S in the Riemannian 
metric defined by the variogram. Specifically, we replace local Euclidean 
distance between points si,S2 £ S hy the square root of the variogram 

(14) Var(Zi(si)-Zi(s2))'/', 

where Zi{s) is the first (say) component of Z{s). Note that the LKC of S 
depends on the local spatial correlation structure of the underlying UGRF, 
although we suppress this dependence in the notation Ci{S). The corre- 
sponding EC density in (13) is then calculated as before, but assuming the 
component UGRFs are isotropic. In other words, the information about the 
nonisotropy of the random field is transferred from the EC density to the 
LKC. 

An explicit expression for the LKC can be found in Taylor and Adler 
[21], but it requires a solid grasp of differential geometry which is beyond 
the scope of this paper. Some idea of how the LKC can be calculated is 
as follows. Suppose S can be embedded by a smooth transformation into 
a set S C M'^ in a higher N > N dimensional Euclidean space so that, in 
this space, local Euclidean distance is the square root of the variogram (14). 
Then 

The Nash Embedding Theorem guarantees that such a finite exists; more- 
over, it is bounded by A < A(A + l)/2 + A. However, in Taylor and Adler 
[21] the actual derivation of the expected EC (2) in the nonisotropic case 
proceeds more naturally from first principles, namely, the variogram metric 

(14) , with isotropic random fields in Euclidean space following as a special 
case. 

Fortunately there is a very simple expression for the A-dimensional LKC, 
which usually makes the largest contribution to the expected EC (2): 

(15) Cn{S)= f det(Var(Zi(s)))i/2^s. 

Js 

Similarly, there is a simple formula for the (A — l)-dimensional LKC: 

CN^iiS) = i / det(Var(ZiT(s)))i/2 ds, 
JdS 

where dS is the boundary of S and Zi(s) is the derivative of Zi tangential 
to the boundary. The lower dimensional LKCs do not have a simple formula, 
except for Co{S) = ^p{S). 

Of particular interest for this paper, it can be shown that in the cases we 
consider, 

C^{SxUd)=fl^{SxUd), 
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where S x Ud is considered as a subset of M + . This justifies carrying on 
with the expected EC (2) as if T{s,u) were jointly isotropic in {s,u). 

FinaUy, Appendix A. 5 gives a neat apphcation of the LKC version of 
the expected EC (13) to derive the scale space EC densities from the 
Gaussian scale space EC densities. 

2.8. The tube method, a direct approach to the P-value of the maximum. 
The tube method, developed by Knowles and Siegmund [10], Johansen and 
Johnstone [9], Sun [15] and Sun and Loader [16], is a direct approach to the 
P-value of the maximum of a random field. 

The approach starts by writing the random field as a Karhunen-Loeve ex- 
pansion, then finding the P-value of the maximum of the first m terms, then 
letting m — > oo. The first m terms can be written as a linear combination 
of m spatial basis functions with m independent N(0, 1) random variables. 
Conditioning on the sum of squares of these Gaussian random variables, the 
P-value of the maximum of the first m terms comes down to the probability 
content of a subset of the m-dimensional unit sphere. This subset turns out 
to be a tube with a radius that depends on the threshold of the random 
field. The P-value calculation is thus reduced to a geometrical calculation: 
the ratio of the measure of the tube to the measure of the sphere. Tak- 
ing expectations over the sum of squares of the Gaussian random variables 
completes the calculation. The final step is to let m — > oo, giving a series 
expansion for the P-value similar to (2). 

Sun [15] gives a general two-term series expansion for the P-value of the 
maximum of a general Gaussian random field that does not have to be 
homogeneous nor have a finite Karhunen-Loeve expansion, and an upper 
bound for an [N — l)-term expansion. In an unpublished manuscript by 
Sun (referenced in Adler [2]), it was generalized to include more general 
boundary cases. Sun, Loader and McCormick [17] developed simultaneous 
confidence regions for response curves in generalized linear models, using 
the tube formula, and generalized inverse Edgeworth expansions, which they 
developed using the Skorohod construction. 

As we can see from this method, there is no need to assume isotropy, 
but the tube method cannot be applied to non-Gaussian fields such as t, 
F or the multivariate random fields that interest us here. Moreover, in the 
case of Gaussian random fields, Takemura and Kuriki [19] show that the 
tube method's series expansion agrees with the expected EC (2) out to the 
number of terms in the expected EC. Thus, when they overlap, the two 
methods give essentially the same result, but the EC method appears to be 
easier to work with, and extendable to the multivariate random fields that 
concern us here. 
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3. Random fields of multivariate test statistics. It is not hard to see how 
we can apply this same union-intersection principle to other types of random 
fields. In Section 3.1 we shall get the EC densities for a Hotelling's field 
from the t field; in Section 3.2 we shall get a P-value approximation (but 
not quite the EC density) for a Roy's maximum root field from the F field; 
in Section 3.3 we shall get similar results for a maximum canonical corre- 
lation field from the correlation field. In fact, all the results so far known 
can be obtained from the correlation field, which makes the computer imple- 
mentation particularly simple. This has been done in the stat_threshold 
function of the FMRISTAT package for the statistical analysis of fMRI data 
(http : //www.math.mcgill . ca/keith/fmristat). 

Further generalizations in Section 3.4 to multilinear forms are obvious. 
In fact, there is a strong link between this paper and Kuriki and Takemura 
[11]. The latter paper is concerned with obtaining P- values for maxima of 
multilinear forms, linear combinations of a multidimensional array (rather 
than just a matrix) of Gaussian random variables. Their interest is only in 
the P- value of the multilinear form itself, not a random field of multilinear 
forms, and their method is based on the tube method described in Section 
2.8. 

The same idea can easily be extended to random fields of statistics 
(Lin and Lindsay [12] and Takemura and Kuriki [18]) by simply replacing 
the sphere Ud with a cone (a subset of the sphere). This will be developed 
in a forthcoming paper [23]. Finally, Appendix A. 5 shows how to get the 
scale space field from the Gaussian scale space field. 

3.1. Hotelling' s T^. Hotelling's field is defined as 

T{s) = ijZ{syW{s)-^Z{s), 

where W{s) is an independent dx d Wishart random field with z/ degrees of 
freedom, generated as the sum of squares matrix of u independent copies of 
Z{s). We are interested in finding pY, the i-dimensional EC density of the 
Hotelling's field. Using Roy's union-intersection principle, we write the 
Hotelling's field as the maximum of the square of a (Student's) t field: 

r(s) = maxr(s, u) where T(s, n) = — ^ \ \ . 
u u W{s)u/iy 

The variance of the derivative with respect to u of Z'u is the identity matrix, 
so that Z'u is a UGRF as a function of u (when restricted to the unit sphere 
Ud) as well as s. Hence, T is the square of a (unit) t field with degrees of 
freedom. 

There is a direct link between the EC of T and the EC of T: for t>0, 
ip{seS: T{s) >t} = ^(p{{s, u)eSxUd: f{s, u) > t}, 
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since for each fixed s, the excursion set {u£Ud- T{s, u) > t} is either empty 
or a pair of ehipsoidal caps with EC equal to 2. To see this, replace uhy v = 
W^^'^u, and note that the excursion set is the intersection of an ellipsoid with 
a pair of symmetric half-spaces. In other words, {{s,u) £ S x Ud'- T{s, u) > t} 
has the same topology as two copies of {s G 5 : T{s) > t}. Figure 3 illustrates 
this point for the simplest nontrivial case of = 1 , d=l and 1^ = 00. 

Let /oj be the fc-dimensional EC density of a t field with v degrees of 
freedom. Then 

E{ip{s G S : T{s) > t}) = iE(v3{(s, u)eSxUd: f{s, u) > t}) 

N+d 

(16) = Yl f^kiS X Ud)pJiVi) 

k=0 

N d 
i=0 j=0 

The first step in the above (16) essentially follows from the expected EC 
(2), though, as noted in Section 2.7, the isotropic theory cannot be directly 
applied. Fortunately, as mentioned in Section 2.7, the Lipschitz-Killing cur- 
vatures oi S xUd agree with the intrinsic volumes oi S xUd considered as a 
subset of M'^"'"'^. In what follows, we will apply this same argument without 
further mention. 

The last step in the above (16) follows from a result of integral geometry 
on the intrinsic volumes of products of sets (see Appendix A. 3). Note that 
the factor of ^ disappears because the expected EC of the excursion set of 
a t field squared is twice that of a t field not squared. Equating (16) to the 
expected EC (2) gives the EC density of Hotelling's random field: 

d 

pfit)=j2p,iUd)pj+,iVt). 

j=0 

Appendix A. 2 gives the intrinsic volume of the sphere pj{Ud), and Worsley 
[26] gives the EC density of the t field. We have thus re-derived the same 
result as Cao and Worsley [4], but with far less trouble. 

3.2. Roy^s maximum root. Let V{s) be a Wishart random field with t] 
degrees of freedom and component random fields independently distributed 
as Z{s). Let Ai(s) > • • • > Arf(s) be the roots of the generalized eigenvalue 
equation V {s)u/ r] = W {s)uXi{s) / v , where W{s) is a Wishart random field 
with degrees of freedom as before. Roy's maximum root random field is 
r(s) = Ai(s). But it can also be derived from the union- intersection princi- 
ple: 

T(s) = max T(s, u) where f(s,u) = )'"'/^ _ 

« u'W[s)u/i' 
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Fig. 3. Example of Hotelling^s field in N = 1 dimensions, with d = 2 components, 
and v — (X degrees of freedom, (a) Excursion set of T{s,u) above t = l (horizontal and 
vertical lines on graphs below); (b) "unwrapped" T{s,u), u — {sm9,cos6) ; (c) Hotelling's 
T^, T{s) = maxn r(s, u); (d) observed and expected EC of T and T as a function of 
threshold t, calculated using (3) and (2). Note that the EC of T is twice that of T (e.g., 
att = l, 6 = 2x3;. 
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Hotelling's random field is a special case with V^(s) = Z(s)Z(s)', that is, 
rj = 1. Unfortunately, the EC of the excursion set of T is no longer directly 
related to that of T, as it was in the case of Hotelling's T^, as the following 
lemma shows: 

Lemma 1. The EC of T is twice the alternating sum of the EC of the 
roots: 

d 

if{{s,u)eS xUd:f{s,u)>t} = 2"^{-iy"^^{seS:Xi{s)>t}. 

1=1 

Proof. The EC of the excursion set of T for fixed s is now more com- 
plicated, but we can find it by Morse's theorem (Morse and Cairns [13]), 
which states the following. Suppose / is any smooth function defined on a 
set A which takes its minimum value everywhere on the boundary of A, and 
all turning points where / = have nonzero det(/) and are interior to A. 
Then Morse's theorem states that 

y'(^) = ^{/{a)=o} sign(det(-/(a))). 

All we have to do is choose a suitable Morse function, /. Let A = diag(Ai, . . . , 
Xd), suppressing dependency on s. Then 

(p{u £ Ud : T{s, u) >t} = {p{a £ Ud ■ a'Aa > t}. 

A suitable Morse function for A = {a£Ud- a'Aa > t} is of course f{a) = a'Aa 
itself. The 2d turning points of / in Ud where / = are the vectors itej, 
where Cj has 1 in position i and elsewhere. The {d — 1) x {d — 1) second 
derivative matrix /(ztcj) is the diagonal matrix with elements \j — Aj for 
all j 7^ i, and so sign(det(— /(±ej))) = (—1)*^^. Adding these contributions 
over the turning points where f >t, that is, Aj > t, gives 

d 

^{u E Ud:f{s,u) > t} = 2^(-iril{;,,>,}, 

1=1 

or in other words, 2 if the number of roots greater than t is odd, and 
otherwise. Applying Morse theory to S x Ud in the same way proves the 
result. □ 

Thus in the case of Roy's maximum root, there is no direct connection 
between the EC of T and the EC of T, except in the case of r] = 1, where 
there is only one nonzero root (equal to Hotelling's T^). For ?] > 1, the 
EC of T is smaller than twice that of T, suggesting that half the EC of T 
approximates the P- value of the maximum root maximized over S: 

P(maxr(s)>t) « iE((^{(s,n) G 5 X Ud:f{s,u)>t}). 
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It is difficult to give quantitative bounds on the error in tfie above approxi- 
mation as tlie fields are non-Gaussian. For Gaussian fields, there exist tools 
such as Slepian's inequality (Adler [2]) and for the finite Karhunen-Loeve 
Gaussian case, the tube work of Taylor, Takemura and Adler [22]. In the 
non-Gaussian case, recent techniques described in Taylor, Takemura and 
Adler [22] give a recipe for bounding the error, though we have not applied 
these techniques here. If the degrees of freedom of the denominator were infi- 
nite, then the maximum of the Roy's maximum root field is a maximum of a 
Gaussian field and the results of Taylor, Takemura and Adler [22] described 
in Section 2.2 apply. 

To evaluate the EC density, note that T is simply an F field with rj, v 
degrees of freedom and fc-dimensional EC density (F^. Then 

E(^{(s, n) G S X i7rf : T(s, u) > t}) = ^ ^^^S)pf{t), 

i=0 

where 

j=Q 

Note that pf {t) is not the EC density of Roy's maximum root; rather, it is 
twice the alternating sum of the EC densities of all the roots (see Figure 4). 
But for high thresholds, the other roots are much smaller than the maximum, 
so their EC is close to zero. For this reason, we can use half p^{t) in the 
approximate P- value (2). 

It is interesting to note that the Roy's maximum root field is not always 
smooth. If the number of dimensions N >2, then it can contain (with posi- 
tive probability) nonsmooth local minima or "cusps" where the two largest 
roots are equal. Figure 5 shows an example with N = 2, d = 2, ij = 6 and 

= oo. The reason is that for equality of the two roots, the two diagonal 
elements vu and V22 of V must be equal, and the off-diagonal element V12 
must be zero. These two constraints are satisfied on two zero contour lines 
of I'll ~ V22 and Vi2- The two lines intersect in points (with positive proba- 
bility) where both constraints are satisfied, and at these points the roots are 
equal. This cannot happen in = 1 dimensions, almost surely, so we do not 
see it in Figure 4. In general, equal roots, and hence cusps, will occur for any 
d>2 whenever N >2. This lack of smoothness rules out the possibility of 
using (8) to find the EC density of the Roy's maximum root field. This does 
not mean a simple expression might not exist; the conjunction random field, 
defined as the minimum of independent smooth random fields (Worsley and 
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Fig. 4. Example of Roy^s maximum root field in N = 1 dimensions, with d = 2 com- 
ponents, and 77 = 2, V = 00 degrees of freedom, (a) Excursion set of T{s,u) above t = l 
(horizontal and vertical lines on graphs below); (b) ''unwrapped'^ T{s,u), u = {sm6,cos6); 
(c) Roy's maximum root T{s) = max.uT{s,u) and the minimum root miiiu T'(s, u); (d) 06- 
served and expected EC of T and observed EC of the maximum and minimum roots, as a 
function of threshold t, calculated using (3) and (2). Note that the EC ofT is twice the EC 
of the maximum root minus the minimum root [e.g. at t = 1, 4 = 2 x (3 — 1)]. Note that 
at high thresholds, the EC of the minimum root is negligible, so the EC of the maximum 
root is well approximated by half the EC of T. 
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1^ > 



Fig. 5. Example of cusps, or nonsmooth local minima, in a Roy's maximum root field, 
with N = 2 dimensions, d = 2 components, and rj = 6, u — oo degrees of freedom. The 
vertical axis is the random field value (the lower lattice is zero). The upper surface is the 
maximum root, the lower surface is the minimum root. Two cusps are visible in the center 
of the picture (above arrow heads) where the upper and lower surfaces almost touch. Cusps 
will occur whenever N >2 with positive probability. 

Friston [30]) also has cusps (at local maxima), yet despite this, a simple 
closed- form expression can be found for its EC density without using (8). 

3.3. Maximum canonical correlation. Let X{r), r G i? C M^^, and Y{s), 
s G 5 C be matrices of UGRFs with c and d columns, respectively, and 
the same number n of rows. Let u £ Uc and v £ Ud- Define the maximum 
canonical correlation random field as 



n 



r,s] = maxTfr, s, n, v) 



where T(r, s, u, v) 



u'X{r)'Y{s)v 



{u'X{r)'X{r)uv'Y{syY{s)vyl'^ ' 

Note that T is the maximum of the canonical correlations between X and 
y, defined as the singular values of {X' X)~^/'^X'Y{Y'Y)~^/'^ . Once again 
there is no direct connection between the EC of T and the EC of T. Using 
the same approach as in Section 3.2, it can be shown that, for positive 
thresholds, the EC of {u, ?; : T(r, s, u, f ) >i} is 2 if the number of canonical 
correlations greater than t is odd, and otherwise. If c = 1 or d = 1, there is 
only one nonzero canonical correlation, so the EC density of T is twice that 
of T. Otherwise, the EC of T is smaller than twice that of T, suggesting 
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that half the EC of T approximates the P-value of the maximum canonical 
correlation maximized over S. This approximation is 

Pf maxr(r,s) > t 

\ r,s 

(17) 

M N c d 

^lY.^'^(^)T.^'^(S)T.^'k{Uc)Y.^'l(Ud)p?,k,J+lit), 

i=0 j=Q k=0 1=0 

where pfj is the EC density of the (cross) correlation random field T for 
fixed u,v given in Cao and Worsley [5]. If M = A^, R = S, and the search 
region in i? x 5" is confined to the "diagonal" where r = s, then this is called 
the homologous correlation random field (Cao and Worsley [5]), which we 
can generalize in the same way. 

We now show that all previous results can be obtained from (17) by setting 
i? to a single point (M = 0), which eliminates the need for special cases and 
makes it much simpler to program. For this reason, and to keep this paper 
self-contained, an explicit expression for pfj is given in Appendix A. 4. First 
of all, it is well known that for fixed r, s the square of the maximum canonical 
correlation is a monotonic function of Roy's maximum root with d as before 
and 

V = Y'X{X'X)-'^X'Y, r] = c, 

W = Y'Y -V, v = n-c, 

and thus, Hotelling's (c = 1) and the F statistic (d = 1) are special cases. 
What is not obvious is that their distributions as random fields are the same, 
since (17) was derived under the assumption of a random X and Y , whereas 
Roy's maximum root (and the others) only require a random Y . However, it 
is easy to see that, conditional on X^V and W have the appropriate Wishart 
distributions with parameters that do not depend on X, hence, they also 
have the appropriate distributions marginal on X. Moreover, since i? is a 
point (M = 0), then X is the same for all s G 5", and so the above remarks 
apply to the random fields, not just the random values at a single s. 

3.4. Maximum multilinear form. Let Uj be a unit dj-vector, j = 1, . . . ,D, 
and u = ui <Si ■ ■ ■ <^ U£i be their Kronecker product, that is, a vector of all 
components of the form ciC2 ■ ■ - cd where Cj is a component of Uj. Let Z{s) 
be a vector of i.i.d. UGRFs of length equal to that of u. The maximum 
multilinear form random field can be defined as 

T(s) = max T{s,ui, . . . ,U£)) where T'(s, ui, md) = Z(s)'ii. 

ui,...,uo 

Note that if Z) = 1, then T(s)^ is a field with di degrees of freedom; if 
D = 2, then T{sf is the maximum root of a di x di Wishart random field 
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with d2 degrees of freedom, the same as v times Roy's maximum root as 
— > oo. The EC density is straightforward: 

PN{t)=(llY. t'k, {Ud, ) p^+fc (t) where k = ^ fe, . 

\j=lfc^=0 / j=l 

Noting that T is unchanged if any pair of uj's are multiphed by —1 suggests 
that 

F maxr(s) > t « 2~^''~^^ ^ 

In the nonrandom field case, = 0, this appears to be the same as that given 
by Kuriki and Takemura [11]; numerical comparisons of the coefficients of 
powers of t are identical in every case considered so far. 



4. Application. 

4.1. Estimating Lipschitz-Killing curvature. Since real data is usually 
nonisotropic, the first step is to estimate the LKC. A method for doing this 
was developed in Worsley, Andermann, Koulis, MacDonald and Evans [29] 
and Taylor and Worsley [24], but, for completeness, we briefly describe it 
here as it pertains to the multivariate linear model (1). 

Let Y{s) = {Yi{s), . . . ,Yn{s))' be the n x d matrix of all the deformations 
data at a point s, and let X = (xi, . . . jXn)' be the n x p design matrix of 
the model (1), assumed of full rank. The least-squares residuals are 

R{s) = Y{s) - X{X'xy^X'Y{s). 

Let Rj{s) be the jth column of R{s), j = 1, . . . ,d. The corresponding nor- 
malized residuals are 

Q,is) = R,{s)/\\R,{s)\\. 

Let Cfc be the A^- vector of zeros with kth component equal to the lattice step 
size along axis k, k = 1, . . . , N . Let 

Dj{s) = {Qj{s + ei) - Qj{s),.. . + cat) - Qj{s)) 

be the n x N matrix of differences of Qj{s) in all lattice directions. Then 
Dj{s)'Dj{s) is an estimator of Var(Zi(s))A, where A is the product of the 
lattice step sizes. Summing over all such lattice points in S and averaging 
over components j, our estimator of the A^-dimensional LKC of S is 

(18) Cj,{S) = Y^Y.^et{D,{s)'D,{s)fl^/d. 
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Worsley et al. [29] show that this is consistent and unbiased in the hmit as 
the mesh size approaches zero. Note that this estimator is not invariant to 
a hnear transformation of the components of Yi{s), and indeed, Rj{s) could 
be replaced by any linear combination of the columns of R{s) and the LKC 
estimator (18) would still be unbiased in the limit. 

Estimating the lower dimensional LKCs is more complicated. Taylor and 
Worsley [24] describe a method that involves filling S with a simplicial mesh 
(a tetrahedral mesh in = 3 dimensions), a nontrivial problem. The coor- 
dinates s € of the mesh are then replaced by the normalized residuals 
Qj{s) £ M". The intrinsic volumes of all the components of the simplicial 
complex (points, edges, triangles, tetrahedra, . . .) are calculated then added 
together in an inclusion-exclusion type formula to give estimates of the 
LKCs of the union of the simplices, and hence, of S itself. Again, this esti- 
mator is consistent and unbiased as the mesh size approaches zero. 

In practice, this estimator produces P-value approximations (13) that are 
a little different from simply assuming that S is an A^-dimensional ball of 
volume Cn{S). The reason is that it is usually the iV-dimensional term that 
makes the largest contribution to the P- value approximation (2). Moreover, 
we have already estimated C]\f{S) quite easily by (18) without filling S 
with a simplicial mesh. In = 3 dimensions, the radius of this ball is r = 
(3/^3 (5) /(47r))"^/^ so that, from (4), the lower-dimensional LKCs can be 
estimated by 

£o{S) = l, A(S)=4r, C2{S) = 2TTr\ 

This short-cut usually results in a slightly liberal P-value approximation 
since a ball has the lowest surface area (and other lower dimensional intrinsic 
volumes) for a given volume. 

4.2. The nonmissile trauma data. As an illustration of the methods, we 
apply the P-value approximations for Roy's maximum root and maximum 
canonical correlation to the data on nonmissile trauma subjects (Tomaiuolo 
et al. [25]) that was analyzed in a similar way in Worsley, Taylor, Tomaiuolo 
and Lerch [31]. The subjects were 17 patients with nonmissile brain trauma 
who were in a coma for 3-14 days. MRI images were taken after the trauma, 
and the multivariate data were the d = 3 component vector deformations 
needed to warp the MRI images to an atlas standard (Collins et al. [6]) 
sampled on a 2 mm voxel lattice. The same data were also collected on a 
group of 19 age and sex matched controls, to give a sample size of n = 36. 
The p = 2 regressors were binary indicators for these two groups of subjects. 
Although a comparison of each trauma case with the control group might be 
useful in clinical applications, we follow Tomaiuolo et al. [25] and pool the 
trauma cases together and compare them with the control group (Tomaiuolo 
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et al. [25] did a similar two-group comparison of white matter density, a 
univariate random field). 

Damage is expected in white matter areas, so the search region S was 
defined as the voxels where smoothed average control subject white matter 
density exceeded 5% (see Figure 1). The LKCs were estimated as in Sec- 
tion 4.1 by those of a ball with volume C3{S) = 2571 (radius r = 8.5). This 
choice of search region is somewhat arbitrary, but a similar search region was 
used in [25] for detecting group changes in white matter density. If it was 
felt that damage was restricted to a smaller region of higher white matter 
density, then the LKCs would be smaller, resulting in lower P-values, lower 
test statistic thresholds and greater power at detecting changes. However, 
this must be offset against the possibility that the smaller search region may 
have excluded regions where change really took place. 

4.2.1. Hotelling^ s T^. The first analysis was to look for brain damage by 
comparing the deformations of the 17 trauma patients with the 19 controls. 
We are looking at a single contrast, the difference between trauma and 
controls, so = 1 and the residual degrees of freedom is = 34. In this case 
Roy's maximum root is Hotelling's T^. The P = 0.05 threshold, found by 
equating (13) to 0.05 and solving for t, was t = 54.0. The thresholded data, 
together with the estimated contrast (mean trauma — control deformations) 
is shown in Figure 1(a). A large region near the corpus callosum seems to be 
damaged. The nature of the damage, judged by the direction of the arrows, is 
away from the center [see Figure 1(b)]. This can be interpreted as expansion 
of the ventricles, or more likely, atrophy of the surrounding white matter, 
which causes the ventricle/ white matter boundary to move outward. 

4.2.2. Roy^s maximum root. We might also be interested in functional 
anatomical connectivity: are there any regions of the brain whose shape (as 
measured by the deformations) is correlated with shape at a reference point? 
In other words, we add three extra regressors to the linear model whose val- 
ues are the deformations at a pre-selected reference point, so now p = 5 
and 1/ = 31. The test statistic is now the maximum canonical correlation, or 
equivalently, the Roy's maximum root for these 77 = 3 extra regressors. We 
chose as the reference the point with maximum Hotelling's for damage, 
marked by axis lines in Figure 1. Figure 1(c) shows the resulting maximum 
canonical correlation field above the P = 0.05 threshold of 0.746 for the 
combined trauma and control data sets removing separate means for both 
groups. Obviously there is strong correlation with points near the reference, 
due to the smoothness of the data. The main feature is the strong corre- 
lation with contralateral points, indicating that brain anatomy tends to be 
symmetric. 



RANDOM FIELDS OF MULTIVARIATE TESTS 



23 



A more interesting question is whether the correlations observed in the 
control subjects are modified by the trauma (Friston et al. [7]). In other 
words, is there any evidence for an interaction between group and reference 
vector deformations? To do this, we simply add another three regressors to 
the linear model whose values are the reference vector deformations for the 
trauma patients, and the negative of the reference vector deformations for 
the control subjects, to give p = 8 and = 28. The resulting Roy's maximum 
root field for testing for these r] = 3 extra regressors, thresholded at 30.3 
(P = 0.05) is shown in Figure 1(d). Apart from changes in the neighborhood 
of the reference point, there is some evidence of a change in correlation at a 
location in the contralateral side, slightly anterior. Looking at the maximum 
canonical correlations in the two groups separately, we find that correlation 
has increased at this location from 0.719 to 0.927, perhaps indicating that 
the damage is strongly bilateral. 

4.2.3. Maximum canonical correlation. If we chose to search over all ref- 
erence points as well as all target points, this would lead to a maximum 
canonical correlation field with X = Y. Note that since the correlation be- 
tween X(r) and Y{s) is the same as that between X{s) and Y{r) (since 
X = Y), then the P- value should be halved. Note also that the reference 
and target points must be sufficiently well separated to avoid detecting 
high correlation due to spatial smoothness. The parameters in our case are 
M = N = c = d = 3, and n is effectively 36 — 2 = 34 after removing the sepa- 
rate group means. The threshold for the maximum correlation random field 
at P = 0.05 is 0.962 from (17) and Section A. 4. Computing all correlations 
is obviously very expensive, but aside from this, the correlation threshold of 
0.962 is so high that the search over all possible correlations is unlikely to 
reveal much beyond the obvious symmetry reported above. 

4.2.4. Conclusion. In conclusion, our analysis shows that damage ap- 
pears to be confined to central brain regions around the ventricles, and not, 
as one might expect, to regions near the brain surface where the brain might 
have impacted the skull. Similar conclusions were reported by Tomaiuolo et 
al. [25] in an analysis of white matter density. 



A.l. Intrinsic volumes. The intrinsic volume ^i{A) of a set yl C M'^ with 
nonempty interior and bounded by a smooth hypersurface is defined as 
follows (Worsley [27]): /Ud(^) is the Lebesgue measure of A, and for j = 
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where C is the curvature matrix of dA, the boundary of A [the {d — 1) x 
{d — 1) second derivative matrix of the inside distance between dA and its 
tangent hyper-plane], detrfc(C) is the sum of the determinants of all k x k 
principal minors of C, and = 27r''/^ /T{k/2) is the Lebes ffue measure of 
the unit {k — l)-sphere in M'^. Note that /Uo(^) = ^{A) by the Gauss-Bonnet 
theorem. 

A.2. Unit sphere. For A = Ud, the unit sphere in W^, fJ,d{Ud) = 0. Now 
C = I(d~i)x{d~i) on the outside of Ud and C = —I(d-i)x{d~i) on the inside. 
Since detvd^i^j{l(d-i)x{d-i)) =2C^J^), then 




_ 2J+V^'/^r((d + l)/2) 

j!r((d + i-i)/2) 

if d — 1 — j is even, and zero otherwise, j = 0, . . . , d — 1. 
A.3. Products. We now show that 

k 
1=0 

A very useful result of Hadwiger [8] states that if ^{S) is a set functional 
that is invariant under rotations and translations, and has the additivity 
property 

(19) ifiSi u ^2) = ^{Si) + ipiS2) - ^{Si n 52), 

then ip{S) must be a linear combination of intrinsic volumes of S. Fixing 
A, we can see that ^{B) = /^^(A x B) has these properties. Fixing B and 
repeating the exercise, we conclude that 

k k 

fj.k{A X B)=J2^Cijfj.i{A)tij{B) 

i=Oj=0 

for some constants Cij. We now determine the constants by judicious choice 
of A and B. First, increasing the size of A by a factor 7 increases its k- 
dimensional intrinsic volume by j^: ^^{'yA) =j'^^k{A). Replacing A,B by 
jA,jB in the above, and noting that x -fB = y{A x B), we conclude 
that the only nonzero constants occur when i + j = k. Next, let A C M*, 
B C M'^"*, both with nonzero Lebesgue measure. Then Ax B cR'' has 
Lebesgue measure ^i{A) ^k~i{B) . Note that ^Xj^A) = for j > i and HjiB) = 
for j > /c — i, so the right-hand side reduces to Ci^k-itJ'i{A)^k-i{B). Since 
Lebesgue measure and intrinsic volume coincide in these cases, we conclude 
that Ci^k-i = 1- This completes the proof. 
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A. 4. EC density of the correlation random field. For completeness, we 
reproduce the EC density pfj{t) of the correlation random field given in 
Cao and Worsley [5]. /OQo(i) is just the upper tail probability of the Beta 
distribution with parameters 1/2, (n — l)/2. Let h = i + j. For i > 0, j > 
and n> h, 

\_{h-l)/2\ 

X ^ |-_-|^^fc^/i-l-2fc|--|^ _ ^2~|(n-l-/i)/2+fc 
fc=0 

/=0m=0^ ^ ^ / \ ^ 

X {l\m\{k — I — m)l{n — 1 — h + I + m + k)] 

X k-l + m)l{j - k-m + /)!)~^, 

where [-J rounds down to the nearest integer, terms with negative factorials 
are ignored and pfj{t) = pfi{t). The summations have been arranged for easy 
numerical evaluation. 



A. 5. Scale space. The Gaussian scale space random field is obtained by 
smoothing white noise with an isotropic spatial filter over a range of filter 
widths or scales, and adding the scale to the location parameters of the field 
(Siegmund and Worsley [14]). In essence, it is a continuous wavelet transform 
that is designed to be powerful at detecting localized signal of an unknown 
spatial scale as well as location. Let dB{s) be Gaussian noise on M'^ based 
on Lebesgue measure and let f{s) be a filter, normalized so that //^ = 1, 
and scaled so that / //' = InxN- The Gaussian scale space random field 
with filter / is defined as 

(20) T{s,w)=w~^l^ I f{{s-t)/w)dB{t). 

Note that T{s,w) ~ N(0,1) and Yar{dT/ds) = w'^^InxN at each point 
s,w. Siegmund and Worsley [14] and Worsley [28] show that for searching 
over a range of scales w E [t(;i,i(;2], 

N 

(21) K{ip{s,w€Sx[wl,W2]■.T{s,w)>t})=Y,^^^iS)pfit), 

i=0 

where the Gaussian scale space EC density is 
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(22) 

i j^^ (1 - 2j)(47ryi!(. - 2j)!^^+i-2^-^ ^ 

[we define w"^ /i as log(it;) when z = 0] . The parameter k measures the variance 
of the derivative in the log scale direction: 

^{s'f + {N/2)ffds, 

with K = N/2 for the Gaussian filter and k = (N + 4) /2 for the Marr filter. 

The scale space random field T{s,w) is nonisotropic in {s,w), so the scale 
space result can be set in terms of the Lipschitz-Killing curvature of 5 x 
ti)2], as in (13): 

Af+l 

(23) K{ip{s,w£ S X [wi,W2]:T{s,w) >t})=Y^ Ci{S x [wi,W2])pf {t). 

Equating the two expressions (21) and (23) for the expected EC implies that 
£0(5' X [t«i,t(;2]) = Ato(5') and for z > 1, 

Ci{S X ['Wi,W2]) 



(24) 



j=0 



X 



i + 2j - 1 

( i-2j)/2(_ip(^ + 2j-l)! 
(l-2j)(4vr).j!(.-l)! 



While the above expression (23) is no more compact than (21), it allows 
us to immediately get an expression for the EC density of the scale space 
random field with d degrees of freedom, defined as the sum of squares of 
d independent copies of the Gaussian scale space field (20). The maximum 
multilinear form random field can be similarly defined. Having expressed 
(21) in terms of Lipschitz-Killing curvatures in (23), it immediately follows 
that the EC density of the scale space field is the same as that of the 
Gaussian scale space field (22), but replacing Gaussian EC densities pf{t) 
with EC densities. We have thus derived the main result of Worsley [28] 
with very little effort. Indeed, we have much more than this result because 
we can compute the EC densities of any scale space random field for which 
we can compute the EC densities in the isotropic, fixed scale, case. 
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